Phase separation at all interaction strengths in the t-J model 
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We investigate the phase diagram of the two-dimensional t-J model using a recently developed 
Green's Function Monte Carlo method for lattice fermions. We use the technique to calculate exact 
ground-state energies of the model on large lattices. In contrast to many previous studies, we find 
the model phase separates for all values of J/t. In particular, it is unstable at the hole dopings and 
interaction strengths at which the model was thought to describe the cuprate superconductors. 



PACS numbers: 71.10.+X, 71.45. Gm, 74.20.-z 

Models idealizing the environment of interacting con- 
stituents are often used to understand collective behav- 
ior. In many cases the problem of providing a solution 
to a model, even though in the model one has abstracted 
the minimal complexity only, is overwhelmingly difficult. 
For example, the antiferromagnetic Heisenberg model, 



Hueis — J ^ ' Si • Sj , 



(i) 



on an infinite square lattice, where the Si are spin-i 
operators, was thought to describe the undoped parent 
compounds of the superconducting cuprates. Despite the 
simplicity of the model, it has no exact solution. How- 
ever, since the model has only spin degrees of freedom, 
numerical solution to the model with satisfactory accu- 
racy has become possible. The physical picture emerging 
from the Heisenberg model and from extensions obtained 
by adding small correction terms to it (which can be un- 
derstood perturbatively) , is in good agreement with the 
experimental results of the undoped cuprates JjJ . 

The next step is to examine the effect of introducing 
holes in a minimal way into a Heisenberg antiferromag- 
net. The t-J model is perhaps the simplest abstraction 
to describe the environment experienced by holes in the 
limit of strong on-site Coulomb repulsion. The i-JHamil- 
tonian, on a square lattice, is written in the subspace with 
no doubly occupied sites as 



■J (7 



J^SvSi-^-), (2) 



where c ia creates an electron of spin a on site i, and 
rii = J2a c la c ic- One can hope, by adding correction 
terms to such a model, some day to understand the doped 
materials. It is clear, however, that no progress can be 
made without the capacity to understand the simple first. 

Even though there are numerous previous studies j^] 
of the t-J model, its phase diagram is unclear. More 
importantly, because of the inherent problem in draw- 
ing unique conclusions from small size numerical calcula- 
tions, the ambiguity and the confusion of the community 
due to such studies grows. One of the main discrepancies 



concerns the phase-separation boundary of the model. 
Emery, Kivelson, and Lin (EKL) || used a combination 
of analytic and numerical calculations to argue for the 
existence of phase separation at all interaction strengths. 
Subsequently many other groups examined larger 

systems numerically and found the model phase sepa- 
rates only for J > t, i.e., outside the physical region of 
the model, and the results of EKL were strongly ques- 
tioned. Most of these studies used a vanishing inverse 
compressibility as the criterion for the onset of phase 
separation Q-^J. The compressibility, however, is not the 
proper observable to find the phase-separation boundary 
in the two-dimensional t-J model, where the transition is 
first order. It is true that the compressibility diverges in 
the region of phase separation, but it jumps discontinu- 
ously at the boundary with the uniform phase. Numer- 
ically, this discontinuity is difficult to see in even large 
finite systems due to the surface energy of the two coex- 
isting phases. The surface raises the energy of a phase- 
separated system, and we find the inverse compressibility 
remains positive even where the system phase separates. 

In this paper we calculate the phase-separation bound- 
ary of the t-J model using the Maxwell construction 
which suffers very little from finite-size effects. We 
present results on the ground states of significantly larger 
size systems than could be studied previously by using a 
new powerful numerical technique. We find a phase dia- 
gram for the t-J model exhibiting phase separation at all 
interaction strengths and in particular in the physical re- 
gion of the model. This result confirms the conjecture of 
EKL and contradicts the more accepted phase-separation 
boundary. We propose a phase diagram of the model 
consistent with all the available reliable results and dis- 
cuss the consequencies of this phase diagram for improved 
models of the copper oxides. 

For sufficiently large J, the t-J model phase separates 
completely due to the attractive nature of the interaction 
term in (^). The two phases have electron densities per 
site of one and zero, respectively. The energy per site in 
the electron region is en — —1.16934 J as determined by 
calculations on the Heisenberg model 0J^]. 

As J/t is reduced, s-wavc electron pairs evaporate from 
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the high-electron density region for J < J c ~ 3.4367i 
P, ^p|JTl] ] . Isolated larger clusters are never stabilized. In 
this range, the two phase-separated regions contain all 
electrons (no holes) and some electrons (some holes). 

To determine the phase-separation boundary for J < 
J c , we use a recently developed Green's Function Monte 
Carlo method to calculate the ground-state energy of 
the model on large finite lattices pp|-p^|. Prior to this 
work, the ground states of correlated fermions on large 
lattices had been calculated exactly only in one dimen- 
sion jl2j or in the limits of small numbers of holes [ p~3|j3~4| 
or electrons jlCj . Starting from an initial trial state, this 
method projects the state iteratively onto its lowest en- 
ergy eigenstate. In principle, this technique could be used 
to project any trial state that overlaps the ground state, 
but trial states with larger overlaps yield smaller statis- 
tical errors and require less computer time. We use gen- 
eralized singlet Jastrow-Slater and Jastrow-paired wave 
functions as initial trial states and have verified conver- 
gence to the ground state by comparing with exact re- 
sults on small systems and by checking that different trial 
states converge to the same state JT3]. Trial states with 
exotic pairing states or broken time-reversal symmetry 
were never needed to converge to the ground state. The 
method suffers from the negative-sign problem, which 
causes the statistical fluctuations to diverge exponen- 
tially with increasing system size at a fixed density. For- 
tunately, with carefully chosen initial trial and guiding 
functions, we are able to calculate the ground state ener- 
gies of significantly larger systems than are possible with 
exact diagonalization. 

To calculate the phase-separation boundary, we need 
the behavior of the ground-state energy as a function of 
density. For a given finite system size, the shape of the 
Fermi surface changes dramatically as the electron num- 
ber is changed, resulting in seemingly random oscillations 
in the energy. To avoid these shell effects, we vary the 
density by changing the system size, keeping the electron 
number and shape of the Fermi surface constant. 

The ground-state energy at J = 2.5t for 32 electrons 
on a variety of system sizes is shown in Fig. These 
finite systems necessarily constrain the electron density 
to be uniform on the length scales of the system size. We 
fit the discrete data to a polynomial, e(n e ), shown as the 
solid curve, in order to treat the energy as a continuous 
function of density. The dashed line, e ps (n e ), is a lin- 
ear function that intersects the Heisenberg energy, en at 
electron density n e — 1 and intersects e(n e ) tangentially 
at a density labeled n ps . 

It is straight forward to show that the ground state of 
the infinite system at a density n e > n ps cannot be a 
uniform phase, because the energy of the uniform phase, 
e(n e ), is higher than e ps (n e ) at the same density. This 
latter energy corresponds to the energy of a mixture of 
two phases, one at electron density ua = 1 and the other 
at electron density ns — n ps . Therefore the infinite sys- 



tem phase separates into two regions with densities ua 
and riB, and its ground-state energy is given by e ps {n e ), 
the value of the dashed line at the average density of the 
system. This is known as the Maxwell construction || . 
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FIG. 1. The ground-state energy per site at J = 2.5i 
for 32 electrons. For clarity, the energies are shifted by a 
linear factor, — e_f/n e . The circles with error bars show the 
energies calculated on lattices of dimensions 6x6, 7x6, 
28 x 28. A sixth-order polynomial fit to the data is shown 
as the solid line, which is extended to the Heisenberg energy, 
the square at energy zero in this shifted plot. The dashed 
line shows the ground-state energy of the infinite system in the 
phase-separated region. We find the onset of phase separation 
occurs at n ps = 0.296±0.004, while the inverse compressibility 
vanishes at n comp = 0.52 ± 0.10. 

The energy of the infinite system is given by the solid 
line in Fig. [j] for n e < n ps and by the dashed line for 
n e > n ps . This Maxwell construction differs from that 
commonly used since the density of one of the constituent 
phases lies at an extreme limit of the allowed density 
range. It is not possible to add electrons to the Heisen- 
berg solid, which has one electron on every site, so the 
dashed line is not tangent to the fitting curve at n e = 1. 
If the t- J model did allow electron densities n e > 1, then 
the intersection point of the solid and dashed lines would 
be shifted to higher densities where the curves could in- 
tersect tangentially. 

In order to be stable, the energy of the infinite system 
must be concave everywhere. Given the solid line in Fig. 
|l| and the allowed density range of the t-J model, the 
dashed line drawn in the figure is the only line possible to 
make the energy of the infinite system globally concave. 

We never examined systems with densities n e > 0.94, 
so we cannot exclude the reentrance of a homogeneous 
phase in this region. For such a phase to be stabilized, 
the solid curve in Fig. |l| would have to drop back below 
the dashed line in this density range. The new Maxwell 
line would lie slightly below the one drawn and would be 
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tangent to the solid curve at both intersections, but the 
phase-separated region would persist at densities n e <■ 
0.94. We never saw any indication of this possibility at 
any interaction strength. 

In the one-dimensional t-J model, the compressibility 
diverges continuously at the transition point in contrast 
to the discontinuous transition in two dimensions 0,ff6|. 
We have verified that in one dimension, the Maxwell con- 
struction yields the same phase-separation boundary as 
that calculated using the inverse compressibility. 




FIG. 2. The phase-separation boundary. The solid line is 
a cubic fit to the phase boundary calculated from systems 
with 32 to 60 electrons as described in the text. Plotted 
are three sample errors, which increase with decreasing J/t. 
The dashed line shows the boundary calculated using just 
8 electrons. The dotted line is a calculation of the limiting 
behavior of the boundary at small electron density. 

Fig. H shows the phase-separation boundary calculated 
at many interaction strengths. Each point is calculated 
from a Maxwell construction using a fixed number of elec- 
trons, either 60, 52, 50, 42, or 32, on at least four differ- 
ent lattices of size L x L or L x L + l where 7 < L < 28. 
At many interaction strengths, we duplicated the calcu- 
lation for different electron numbers, and always found 
the discrepancy in the critical density to be comparable 
to the statistical errors. We have no reliable results for 
J < O.li, where the phase-separation boundary extends 
beyond our maximum electron density of n e « 0.94, cal- 
culated using 60 electrons on an 8 x 8 lattice. 

At small electron density the phase-separation bound- 
ary is given by J c (n e — >• 0) = J c {n e = 0) + j^tn e , which 
is plotted as the dotted line in Fig. |[ This expression is 
obtained by assuming that the kinetic energy of paired 
electrons varies with density while the binding energy is 
invariant. The calculated boundary approaches this slope 
at low electron density. 

The dashed line in Fig. shows the phase-separation 



boundary calculated using the Maxwell construction with 
8 electrons on 4 x 3, 4x4, and 5x4 lattices. The 
energies of these small systems may be calculated with 
Green's Function Monte Carlo or exact diagonalization. 
The discrepancy between these results and those using 
32 or more electrons is apparent, but both boundaries 
extend to very small J/t. 

The numerical calculation of EKL varied the electron 
density of a 16-site lattice, resulting in both finite-size 
and shell effects absent in the present work ||. However, 
they also found phase separation at all J/t. 

We now use our knowledge about the various regions of 
the phase diagram to derive a proposal for the full phase 
diagram of the model, shown in in Fig. ^. We know with 
high degree of confidence the phase-separation line ob- 
tained in this work, shown by the solid line. Wc also know 
what happens (almost rigorously) at low electron density 

@,[y]|i7m]. In thc limit of 

zero density, electrons form 
s-wave pairs for J > 2i, and these pairs phase separate 
at J > J c ~ 3.4367t. For infinitesimal densities, the elec- 
trons are unstable to higher-angular-momentum pairings 
due to the Kohn-Luttinger effect [ pj| . The strongest in- 
stability, as determined by T- matrix calculations, is p- 
wave at small J/t and d 2 .2_ y 2-wave at intermediate in- 
teraction strengths [ pd)Jl~^ ]. We have drawn the bound- 
aries between these phases as solid lines at small densities 
where the calculations are valid. We believe these phases 
continue to higher densities but cannot trust the precise 
location of the lines obtained by the low density expan- 
sion, so we extend the boundaries as dashed lines. We 
believe the lower boundary of the d-wave phase meets the 
phase-separation boundary at J* ~ 0.27t, the minimum 
interaction strength for which a two-hole d-wave bound 
state is stabilized |l4j . With the present work's evidence 
of phase separation, we interpret the two hole binding as 
indicating phase separation into a state of d-w&ve sym- 
metry. Therefore, we extend the lower boundary of the 
d-wave phase using the expressions from Ref. (lj] shifted 
linearly to extrapolate to J*. 

For J < 0.1 the model has a ferromagnetic instability, 
which is sketched schematically as the dotted line |l9) . At 
low hole dopings, this boundary obeys en = 2tt(1 — n) 2 t, 
so J/erro — 5.37(1 — n) 2 t || . Series expansions indicate 
this phase does not extend beyond n e < 0.7 [M. In 
the unpolarized region, a continuation of the low density 
p-wave pairing phase is compatible with ferromagnetic 
correlations near the ferromagnetic instability. 

This phase diagram will be sensitive to other terms 
that one may add to the t-JHamiltonian. In particular, a 
long-range Coulomb repulsion will suppress macroscopic 
phase separation. The local tendency for phase separa- 
tion, however, could have consequencies for the dynam- 
ics of the more complete Hamiltonian. The competition 
of this tendency with the Coulomb forces might create 
stable clusters with definite sizes, as in nuclear droplets 
(nuclei) , or it might give rise to either the "local" collec- 
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tive modes, such those considered by Emery and Kivelson 
pl[ , or to some superlattice structure, such as those seen 
in neutron diffraction experiments p2| . 
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FIG. 3. Our proposed phase diagram. The boundaries 
drawn as solid lines are accurate, while the others are more 
speculative. The phase-separation curve is the fit to the re- 
sults shown in Fig. ^| The boundaries between phases of 
different pairing symmetry are calculated from expressions 
in Ref. JIT]]. They are accurate at low electron density, where 
they are drawn as solid lines, and are extended as dashed lines 
to the phase-separation boundary. There is a ferromagnetic 
phase at very small J/t and small hole dopings. 

In conclusion, we have shown the two-dimensional t-J 
model phase separates at some range of densities for 
J > O.lt, and we believe the instability extends to all pos- 
itive interaction strengths. The only assumption made is 
that the uniform phase is unpolarized in this region. In 
particular, phase separation occurs in the region of pa- 
rameter space where the model is thought to apply to the 
cuprate superconductors. We verified that phase separa- 
tion extends to electron densities of at least n e rj 0.94 
at all interaction strengths, and we believe it extends to 
n e = 1, the undopcd antifcrromagnet. The main rea- 
son for the discrepancy of these results with previous 
work is that the phase-separation boundary is determined 
far more accurately using the Maxwell construction than 
the inverse compressibility. Finally, we have proposed a 
complete phase diagram of the model including pairing 
symmetries based on all accurate calculations presently 
known to us. 
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